Universal three-body recombination via resonant d-wave interactions 
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For a system of three identical bosons interacting via short-range forces, when two of the atoms 
are about to form a two-body s-wave dimer, the Efimov effect takes place leading to the formation 
of an infinite number of three-body (Efimov) states. The lowest Efimov state crosses the three-body 
break-up threshold when the s-wave two-body scattering length is a ~ — 9.73r v( jw, »"vdw being the 
van der Waals length. This article focuses on a generalized version of this Efimov scenario, where 
two of the atoms are about to form a two-body d-wave dimer, resulting in strong d-wave interactions. 
Bo Gao has argued, in [Phys. Rev. A. 62, 050702(R) (2000)], that for broad resonances the d-wave 
dimer is always formed when a ~ 0.956r v dw- Our results demonstrate that a single universal three- 
body state associated with the d-wave dimer is also formed near the three-body break-up threshold 
at a ~ 1.09r v dw, or alternatively 02 = 0.902r v dw, where 0,2 is the two-body d-wave scattering length. 
Such a universal three-body state is signaled experimentally by an enhancement of the three-body 
recombination rate. The three-body effective potential curves that are crucial for understanding the 
recombination dynamics are also calculated and analyzed. An improved method to calculate the 
couplings, effective potential curves, and recombination rate coefficients is presented. 



I. INTRODUCTION 

In recent years, the investigation of Efimov physics [l[ 
has attracted much interest and has sparked a tremen- 
dous advance in our understanding of its fundamental 
aspects from both experimental |3-[l2| and theoretical 
[13Ml8l | viewpoints. Efimov physics describes a universal 
phenomenon in three-body systems: when the atoms are 
about to form a two-body s-wave dimer, i.e., the s-wave 
scattering length a goes to infinity, and an infinite num- 
ber of three-body Efimov states exists. The ratio between 
the energies of nearby Efimov states is given by a univer- 
sal scaling factor E n+ i/E n = e ~ 27r / s ° , where so depends 
only on the mass ratios, number of resonant interactions 
and quantum statistics of the atoms, e.g., e*l Sa ps 22.7 
for three identical bosons. 

Therefore, a single parameter, called the three-body 
parameter, is needed to determine the whole Efimov spec- 
trum and other scattering observables. One possible def- 
inition of the three-body parameter is a*_ , the value of a 
at which the first resonance in the three-body recombi- 
nation rate K3 appears, i.e., at which the lowest Efimov 
state energy crosses the three-body break-up threshold 
(E = 0). One fundamental assumption underlying all 
Efimov physics is that the three-body parameter encap- 
sulates the short-range details of two- and three-body 
interactions, and thus should not be universal. Surpris- 
ingly, however, ultracold experiments with alkali atoms 
near Fano-Feshbach resonances have observed a 

universal value, namely al ~ — 9.1r v dw m homonuclear 
atomic systems, where r v dw is the van der Waals length. 
The van der Waals length r v dw = i^^bCe) 1 ^ 4 /2 is a 
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characteristic length scale for the van der Waals inter- 
action —Ce/r 6 between two neutral atoms with two- 
body reduced mass /i 2 h. A relevant energy scale, called 
the van der Waals energy, can be defined as i? v dw = 
h 2 / (2/i 2f ,r 2 dw ). 

This newly found universality in the three-body pa- 
rameter was sub seq uently studied in several different the- 
oretical models [ljJUl]. In particular, Refs. [22|, [23| con- 
clude that the universality of the three-body parameter 
comes from a universal barrier in the three-body poten- 
tial curves arising from interatomic distances around 2 
r v dw • This barrier originates from the suppression of the 
probability to find two atoms at distances shorter than 
r vc jw due to the increase of the classical local velocity and 
it is a universal property of van der Waals interactions be- 
tween neutral atoms. Our previous studies [22|, [24[ have 
indicated that the Lennard- Jones potential is an excellent 
model potential to study the universality of three-body 
physics in ultracold atomic gases, and for this reason it 
is also adopted in this analysis. 

The present study, however, employs the Lcnnard- 
Jones potential model to study three-body systems with 
atom pairs close to a d-wave resonance. Interestingly, 
Gao [25[ predicted that d-wave dimers form at a universal 
value of a, a rj 0.956r v dWj for potentials with a van der 
Waals tail. Evidently, since that prediction was obtained 
for a single channel potential model its validity is ex- 
pected to hold only for broad Fano-Feshbach resonances, 
and possibly is limited only to single-channel two-body 
interactions. Nevertheless, one natural question is: will 
there also be a universality of three-body physics near 
that value of scattering length, i.e., when d-wave inter- 
actions are expected to be strong? Here, the three-body 
recombination rate is studied to answer this question, and 
it is seen to exhibit two universal enhancements. Anal- 
ysis of the three-body effective potential curves yields 
an intuitive understanding of these enhancements. One 
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of these enhancements is caused by the simple fact that 
when the two-body d-wave state becomes barely bound 
an additional decay channel is formed, leading to an en- 
hancement in recombination to d-wave dimers. The sec- 
ond enhancement, however, has a different nature. It 
corresponds to the formation of a universal three-body 
state that crosses the collision threshold for three free 
particles. Our analysis indicates that other three-body 
states associated with excited 2-body angular momenta 
are also possible, however, their signatures in recombi- 
nation will be suppressed in the low energy limit due to 
the Wigner threshold laws for recombination [26[. It is 
important to mention that since in experiments using ul- 
tracold atoms s- and d-wave states are coupled by the 
magnetic anisotropic dipole interaction the s-wave scat- 
tering actually diverges when the d-wave dimer becomes 
bound. Although the connection between these three- 
body states associated with d-wave resonant interactions 
and the multichannel physics in ultracold atoms still re- 
mains to be understood more deeply, the present study 
can offer an alternative parameterization in terms of the 
properties of the d-wave scattering length. Moreover, 
since the three-body effective potential curves are an im- 
portant tool for understanding universality in few-body 
systems at ultracold temperatures [22T - |24j ] . an improved 
method to calculate the coupling matrix elements is also 
developed in this paper. 

The remainder of the paper is organized as follows. 
Section II tests the Gao prediction numerically for a 
Lennard-Jones potential, and demonstrates very good 
agreement. The three-body recombination rate is then 
calculated using the Lennard-Jones two-body potential 
with 2 s-wave bound states in the vicinity of the d-wave 
resonance. The improved numerical method is discussed 
in the subsection II. A. Analysis of the recombination rate 
using effective potential curves is made in subsection II. B. 
Section III checks the universality using the Lennard- 
Jones potential with 3 s-wave dimer bound states. 



II. UNIVERSAL TWO-BODY d-WAVE DIMER 

One of the key experimental tools used to study ul- 
tracold atomic gases is the Fano-Fcshbach resonance, 
which can be utilized to magnetically tune the scatter- 
ing length. Although the multichannel character of the 
hyperfme interactions leads to a great deal of complica- 
tion, the single-channel van der Waals interactions has 
been shown to offer a good model for studying a broad 
Fano-Feshbach resonance in ultracold atomic gases [27j . 
In 2000, Bo Gao predicted that for single-channel inter- 
actions with a van der Waals tail, — Cg/r 6 , there is always 
a d-wave dimer (and dimers with higher angular momen- 
tum I = 4j + 2, where j = 1,2,3...) that becomes bound 
at a universal value of the two-body s-wave scattering 
length a = a* = 47r/[r(l/4)] 2 r vdw w 0.956r vdw j^. 
[Note that a* has the same value as the so-called mean 
scattering length a as defined by Gribakin and Flambaum 



[271 . [281 ] , which in turn should not be confused with the 
usual scattering length a.] Another special property of 
the two-body system at a = a* is that the effective range 
r e ff has the value 2[r(l/4)] 2 / (3ir) r vdw rj 2.789r vdw , 
which equals the effective range when s-wave scattering 
length is infinity. In Gao's work [25j . the universality 
of a* is based on an /-independent quantum-defect the- 
ory, where the short-range physics can be described by a 
single short-range parameter, K c , which is itself approx- 
imately /-independent under proper conditions that are 
satisfied by most diatomic systems. In combination with 
solutions for the long-range interactions, K c uniquely de- 
termines both the scattering and the bound state prop- 
erties of diatomic systems. 
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FIG. 1: (Color online) The values of the two-body s-wave 
scattering length a,* at the point where a d-wave (I = 2) dimer 
(black curve with square symbols) and where an i-wave (I = 6) 
dimer just becomes bound (red curve with circular symbols) 
are shown as functions of the number of two-body s-wave 
bound states. 

In this section, the universal value of the two-body s- 
wave scattering length a ; * [where a d-wave (I = 2) or 
i-wave (I = 6) dimer becomes bound] is studied nu- 
merically MS cl function of the number of two-body s- 
wave bound states and using the Lennard-Jones potential 
model for the two-body interaction: 



(i) 



Note that in the present work the parameter A is used 
to adjust the values of a as well as the number of bound 
states. (Some additional details about the Lennard-Jones 
potential are discussed in Appendix [5]) In Fig. [T] the 
black square symbols (red circular symbols) indicates 
the numerical values of the two-body s-wave scattering 
length a* 2 (ag) at the point where a d-wave (i-wave) dimer 
becomes bound. The more s-wave bound states sup- 
ported by the system, the deeper is the potential, im- 
plying different short-range physics. The universal pre- 
diction from Gao's work, a* « 0.956r vd w [H[ (horizontal 
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dashed line in Fig. [IJ, however, suggests that the actual 
form of the potential at short-distances is not important. 
In fact, our numerical results obtained for a ; * agree well 
with the Gao prediction of two-body binding energies 
[within 1% (6%) in the case of 10 s-wave bound states 
for a d-wave (i-wave)] and this agreement tends to im- 
prove as one increases the number of s-wave states. In 
the vicinity of the point where the dimcr is just barely 
bound, i.e., when a < a*, the binding energy can be ex- 
pressed as a linear function of the scattering length , i.e., 



Eij [fr 2 /(2^2b^ dw )] ~ d\ {a\ - a) /r vdw , 



(2) 



where di « 5.8 and d± ss 43, for d-wave and i-wave states, 
respectively, are approximately universal. It is also inter- 
esting to note that near a = a* 2 one expects strong d-wave 
interactions. In order to parametrize d-wave interactions 
in the ultracold regime, scattering properties must also 
be explored for two atoms colliding in the I = 2 channel. 
In this case, due to the 1/r 6 dependence of the Lennard- 
Jones potential, the d-wave elastic scattering phase-shift 
82 (k) is known to have the following low energy expan- 
sion: 



tan ($2 (k) = — 



5 ;„5 



a° 2 k 



(3) 



instead of the usual k 2l+1 dependence found for I < 2. 
Here, A 2 - [71T (5) T (1/2)] / [4T 2 (3) T (11/2)] r^ dw a 
0.160r^ dw is a constant that only depends on the van 
der Walls length, and 02 diverges when a d-wave state is 
just about to be bound (2{|. Therefore, 0,2 is denoted the 
"d-wave scattering length" in this study (which should 
not be confused with a%), and is also used in order to 
characterize our findings in terms of the d-wave interac- 
tions. This parameterization, therefore, allows a better 
comparison to experiments in ultracold quantum gases. 
It is important, however, to emphasize that the com- 
plicated multichannel structure of alkali atoms used in 
ultracold quantum gases can strongly affect the relations 
explored above between the three-body states found for 
d-wave interactions and the single channel parameters 
that characterize the present study. This is particularly 
true near narrow Feshbach resonances 11411. 



A. Hyperspherical approach 

The hyperspherical approach is utilized here to study 
the three-body systems. The version of that approach 
adopted here has been discussed in detail in other recent 
studies [22|, [3(j. To briefly summarize, the three-body 
system is described in hyperspherical coordinates: a hy- 
perradius R that describes the overall size of the three- 
body system, two hyperangles 9 and <fi that describe the 
shape of the three-body system, and the three usual Eu- 
ler angles a, /3, 7 that describe rigid-body rotations. In 
these coordinates, the three-body Schrodinger equation 
is written as, 



/ d 2 A 2 + 15/4 



+ V(R,6,tp)-E 



01 



0. 



2/x 36 \dR 2 R 2 

(4) 

where A 2 = A 2 (ft) [ft = {9, (p, a, /3, 7}] is the "grand 
angular- momentum operator" and fi^b = is the 

three-body reduced mass of three identical atoms with 
mass m [3l|, |32| ■ 




FIG. 2: (Color online) Adiabatic three-body potential curves 
for a ~ 0.977r v dw- The 2 + 1 channels, which corresponding 
to a dimer plus a free atom at very large distance, are labeled 
by a combination of a letter and a number. The letter denotes 
the angular momentum quantum number / of the dimer, and 
the number labels the channels for the same dimer angular 
momentum from low-to-high dimer binding energies. 



III. THREE-BODY RECOMBINATION 

This section focuses on the three-body recombination 
rate near the region where the d-wave dimer is about to 
be bound. In particular, the example of a Lennard-Jones 
potential with two s-wave bound states is elaborated. In 
this numerical example, the d-wave dimer becomes bound 
at about a = 0.995r vd w and the i-wave dimer becomes 
bound at about 1.206r vd \y- Later, the three-body recom- 
bination rate with more two-body s-wave bound states 
will be computed to test the universality of those results. 



To solve this Schrodinger equation, the first step in 
our method is to diagonalize the adiabatic Hamiltonian 
iJ a d {R, ft) and find the adiabatic potentials U V (R) and 
channel functions $„ (R; ft): 



Had (R, ft) *„ (.R; ft) - U v (R) (R; ft) 



(5) 



whose solutions depend parametrically on R. The adi- 
abatic Hamiltonian, containing all angular dependence 
and interactions, is defined as 



Had (-R, fi) 



ft 2 A 2 



15ft 2 



2fi 3b R 2 8n3bR 2 



+ V(R,9,ip) 



(G) 
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For illustration, Fig. [2] shows the three-body adia- 
batic hyperspherical potential curves for the case of a rj 
0.977?Vdw- [Note that in this figure the labels for the 
final recombination channels are also indicated and are 
useful for further analysis.] 

The adiabatic potentials U V (R) and channel functions 
$„(P; 0) obtained by solving Eq. for fixed values of 
P contain all the correlations relevant to this problem. 
For each P, the set of $ v (O; R) is orthogonal, 

[ dn<f> ll (R;n)*$ u (R;n) = 6^, (7) 

and complete 

^ $ T (i?; ft) $ T (i?; IT)* =5(0-0')- (8) 



In practice, calculation of all the channel functions is time 
consuming and impractical. However, numerical studies 
have shown that only a small number of channels are 
needed as a truncated basis-set to expand the total wave 
function, e.g., 



N 



ii E (R,n)=Y^F®(R) (0;P), 



(9) 



where N c is the number of channels adopted. Insertion 
of Eq. © into Eq. (TJ| leads to a set of coupled one- 
dimensional equations: 



2^ 3 fc dR 2 



U v (R) - E 



E 



r 



2P„„ (R) 



dR 



Qvn (R) 



F E , 



(R) = o, 



(10) 



where P and Q are the coupling matrices defined below, 
and v' denotes the v 1 — th independent solution. (Here- 
after, unless otherwise specified, we use an underline to 
denote the matrix form, e.g., P denotes a matrix with 
matrix element P,^.) In the equation above, the non- 
adiabatic coupling matrices are defined as 



P vll (R) 



dfM>„(P;0)* A$ M (P;0), (11) 



potentials are usually more physical than the raw Born- 
Oppenheimer potentials because they include hypcradial 
kinetic energy contributions through the P 2 V term. For 
example, the effective potential gives physical asymptotic 
behaviors of the system at large R with the correct coef- 
ficient of R~ 2 . In some cases, one can even neglect the 
nonadiabatic couplings between different adiabatic chan- 
nels and approximately solve the Schrbdinger equation 
as a single channel problem: 



Qv/j, (R) = J dn$„ (R; 0)* — $ M (P; O) , (12) 

and are responsible for inelastic scattering processes as 
well as the finite width of three-body resonant states. In 
practice, only the 



(13) 



component of Q v ^ is needed to solve the coupled equa- 
tions [Hj]. The relation between Q and P 2 is given by 
^P = —P^_+Q. From the definition of the P and Q ma- 
trices, it is easy to see that the coupling matrices have the 
following properties: P„ M = — P^ and P 2 = P 2 „, which 
leads to P vv — 0, and Q vv = P 2 „. In the hyperspheri- 
cal adiabatic representation, the effective potentials are 
usually defined as 

U v (P) = W vv (P) = U v (P) - (R) , (14) 

and differs for the hyperspherical potential U v (R) by the 
addition of the diagonal correction P 2 „(P). The effective 



2^36 dR 2 



U v (P) - E 



Ft (R) = 0, 



(15) 



called the adiabatic approximation [34J. Although a 
semi-quantitative picture can emerge within this approx- 
imation, the fully coupled system of equation has to be 
solved to obtain a more quantitative description of the 
properties of the three-body system including inelasticity 
in particular. 

Quantitative solutions of Eq. (| 10[) hinge on the ability 
to calculate non-adiabatic couplings accurately. A tradi- 
tional method for calculating the coupling matrices is to 
apply a simple finite difference scheme for the derivative 
of $ M (P;0), i.e., 



or ^ 



(P;0) 



(P + AP; O) - % (P - AP; O) 
2AP 



(16) 

This scheme, however, is only accurate up to the first 
order in AP. Here, we use instead an improved method 
to calculate the 9$ M (P; O) /dR. This method only relies 
on P ad (P, O) and dH ad (P, O) /dR at the desired value 
of P and is shown to be in principle numerically exact 
(see Appendix [Bl for a detailed discussion). 
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B. Three-body recombination rate 

The numerical method described in the previous sub- 
section yields more accurate coupling matrix elements 
and hence improved effective potential curves. These ef- 
fective potentials are not only crucial for intuitively un- 
derstanding the system, but also for quantitative calcu- 
lations of three-body observables like the recombination 
rate. Together with the i?-matrix propagation method 
discussed in Ref. [3(| , the three-body recombination rate 
for J v = + symmetry, where J is the three-body total 
angular momentum and ir the parity, are calculated and 
shown in Fig. [3] The solid black curve with square sym- 
bols shows the total recombination rate as a function of 
the scattering length a, near the d-wave resonance, i.e., 
for values of a near a\ . The total rate exhibits two clear 
enhancements labeled by "A" (oa ~ l-09r v( jw) and "B" 
(ob ~ 0.98r v dw)- The enhancements "A" and "B", indi- 
cated in Fig. [3] by the vertical dashed lines, occur right 
after and right before the d-wave dimer becomes bound 
(see vertical solid line in Fig. \5§ , respectively. Analysis 
of the partial recombination rates helps to elucidate the 
mechanism responsible for these enhancements. The en- 
hancement around the peak "A" appears in all the par- 
tial three-body recombination rates. After the d-wave 
becomes bound and forms a new decay channel, the par- 
tial J<3 rates near the enhancement "B" have a different 
characteristic. The total rate is now dominated only by 
the partial rate into this new d-wave channel [see the inset 
of Fig. [3] . These features suggest that peak "A" might 
be a resonance due to a three-body state that crosses the 
three-body threshold, while peak "B" seems to relate to 
some other threshold behavior related to the formation of 
the new d-wave decay channel. The effective potentials 
near these enhancements, however, can yield insights into 
their origin. 

First, consider the physical origin of the enhancement 
A. Figure. [4] (a) shows the adiabatic hypcrspherical po- 
tential curves at peak A. Only the channels relevant to 
the resonances are shown in this figure. The red dashed- 
line shows the effective hyperspherical potential. We have 
diabatized the potential near an avoided crossing since it 
plays no role, and display it as the black solid curve. The 
black solid curve shows an outer barrier, and a three- 
body short-ranged state could be supported for this po- 
tential. To check whether there is a three-body state 
near the threshold, a calculation of the WKB phase at 
zero scattering energy is carried out: 

0wkb = / =dR, (17) 

Jr. y/-2» 3b W vv (R) 

where R a and Rb are the classical turning points shown 
in Fig. S] (a) . The calculated WKB phase </>wkb is about 
0.517T, which is strong evidence for the existence of a 
three-body state. [Below, we calculate the actual en- 
ergy of this state and describe some of its properties.] 
It is this three-body state associated with the d-wave 




FIG. 3: (Color online) Total and partial J n = + three-body 
recombination rates as functions of the two-body scattering 
length a, shown on a linear scale, near a d- wave resonance 
associated with the d-wave dimer formation. The inset shows 
the same graph with a logarithmic scale for the y-axis. The 
units of Kz are converted to cm 6 /s by using the van der Waals 
length r vdw = 101.0 bohr and mass 132.905429 amu of 133 Cs. 
The solid vertical line indicates the value of the s-wave scat- 
tering length at the point where the d-wave dimer becomes 
bound. The two dashed vertical lines indicate the two values 
of a at which recombination rate is enhanced, denoted A and 
B respectively. 

dimer state crossing the three-body break-up threshold 
that causes the enhancement A in the three-body recom- 
bination rate. The physical picture is that the three-body 
state at threshold can trap the system within the poten- 
tial well for a long time, hence making it more likely that 
the system will decay into deeper channels. Therefore, 
this decay process only depends on the initial channel, 
which is consistent with the fact noted above that near 
the enhancement A all partial rates display the resonance 
effect. 

Next consider the enhancement mechanism responsible 
for peak B. Figure 2] (b) shows the effective potentials 
W VI/ (R) as functions of R for three different scattering 
lengths around the enhancement peak B as solid curves, 
dashed curves, and dash-dotted curves. The correspond- 
ing K3 rates are indicated by the three points (£), (ii), 
and (Hi) in the inset of Fig. 3) Here, the dominant path- 
way at low energy is from the lowest entrance channel 
to the newly formed d-wave channel. For each scatter- 
ing length, the thin curves in Fig. 2] are the entrance 
channels, and as one can readily see, they do not change 
substantially in all three cases. This means that the inci- 
dent channel potential curve depends on the atom-atom 
scattering length only weakly. The thick curves in Fig. @] 
are the most important recombination channel associated 
with the new d-wave channel. It is clear that from point 
(i) to (Hi), the potential barrier decreases appreciably 
and can explain the increase of recombination as it grad- 
ually allows the three-body wave function to approach 
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shorter distances, where non-adiabatic couplings are ex- 
pected to be stronger. 

From Ref . (30j , the scattering matrix element describ- 
ing the recombination from the lowest incident channel 
to the newly formed <i-wave channel can be approximated 
by, 



2m / dRF f (R)W f i{R)Fi(R) (18) 



where Fi {Ff) are the regular solutions of Eq. (|15|) — sin- 
gle channel approximation — for the initial (final) chan- 
nel i (/). As mentioned above, the incident channel is 
dominated by the centrifugal barrier which, in turn, is 
independent of the scattering length. Hence the energy- 
normalized Fi are well approximated by Bcsscl functions: 



Ft (R) « vW?^2 (hR) , 



(19) 



which is, of course, also independent of scattering length. 
The coupling matrix element Wfi (R) is also found, nu- 
merically, to depend on scattering length only weakly for 
the three cases we studied here. Therefore, we assume 
that the product of Wfi (R) Fi (R) does not depend ap- 
preciably on the scattering length. In addition, this prod- 
uct is shown numerically to be significant only at small 
hypcrradius (around 1.3 r v( jw to 2.8 r v aw)- Therefore, 
the scattering length dependence of Sfi is fully deter- 
mined by Ff {R) at short distances. There are two com- 
peting mechanisms controlling the amplitude of Ff (R): 
one is the tunneling through the outer barrier, which 
makes the amplitude become increasingly larger, as the 
barrier becomes increasingly lower from (i) to {Hi); the 
other is the overall amplitude for an energy normalized 
wave function, which makes the amplitude becomes in- 
creasingly smaller, since the threshold energy associated 
to the energy of the d-wave dimer become lower and 
lower, and the out-going wave vector become larger and 
larger. These two competing mechanisms help to form 
peak structure B. The WKB phase calculated for these 
three cases is smaller than 7r/2, which is another evi- 
dence that this enhancement is not due to a three-body 
state resonance. This is also consistent with the fact that 
such an enhancement is not shown in the partial recom- 
bination rates for channels other than the newly formed 
d-wave. 



IV. THREE-BODY STATE ASSOCIATED WITH 
THE d-WAVE DIMER 

While the Efimov states can be viewed as three-body 
states associated with a s-wave dimer near the three- 
body breakup threshold, the new three-body state dis- 
cussed in the previous section is a three-body state that 
can be associated with a d-wave dimer. Because of the ex- 
istence of deeper atom-dimcr thresholds, this three-body 
state is actually a quasi-bound state. Figure [S] shows the 
three-body state energy as a function of the scattering 
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FIG. 4: (a) (Color online) Ww{R) as a function of hyperra- 
dius R at the Kz resonance peak A. (b) (Color online) W{R) VV 
as a function of hyperradius R around the K3 enhancement 
peak B. The inset shows the three points where the potential 
curves correspond to. 



length a (red dots), with the error bars indicating the 
width of the quasi-bound state. The red line represents 
a fitting formula to this energy: 



E, 



h 2 / (2/x 2b r 2 



= d. 



36" 



vdW) 



( a *3b - a ) 
TVdW 



e3b 



(«3& " a Y 



(20) 



1 vdW 



where a^ b w 1.09r v dWi ofab ~ 8.21 and e^b ~ 10.01 are 
fitting parameters. For comparison, the figure also shows 
the d-wave energy as the black curve with solid square 
symbols. The three-body energy crosses the three-body 
breakup threshold at a ~ 1.09r v dWi which corresponds 
to peak A in if 3. As one can see, the three-body state is 
formed before the d-wave state becomes bound and its en- 
ergy is always below the d-wave dimer energy (the black 
curve with square symbols in Fig. [5]). The difference be- 
tween the dimer and trimer energy also increases when 
the scattering length becomes smaller, and the dimer be- 
comes deeper. However, the width of the three-body 
quasi-bound state docs not seem to change appreciably. 
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FIG. 5: (Color online) Energy of the three-body bound state 
associated with a d-wave dimer as a function of scattering 
length a, both in van der Waals units. 



FIG. 7: (Color online) The total three-body recombination 
rates as a function of r^ dw /a|. The units of Kz are the same 
as Fig. [6] 



The Gao analysis in Ref. [25[ , suggests that there is also 
an i-wave dimer that becomes bound at a nearby scat- 
tering length. However, the i-wave dimer does not affect 
these two peaks. The i-wave channel has a very sharp 
avoided crossing with all the other channels. As a conse- 
quence, the partial K3 rate into the i-wave dimer channel 
is negligible compared with the other partial rates. One 
explanation may be that forming a i-wave (/ = 6) dimer 
results in the exchange of a large amount of angular mo- 
mentum between the dimer and the free atom. There- 
fore, although higher partial waves I = 10, 14, 18... might 
also be formed at nearby scattering lengths, they are not 
expected to show any strong features in the three-body 
recombination rates at ultracold energies. 
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LJ with 2 s-wave 
LJ with 3 s-wave 



1.0 1.1 



1.2 1.3 



1.4 



FIG. 6: (Color online) The enhancements for the total 
three-body recombination rates at about a — 0.995r v dw for 
Lennard- Jones potential with 2 and 3 s-wave bound states. 
K3 is convert to cm 6 /s by using van der Waals length r v( jw = 
101.0 bohr and mass 132.905429 amu of 133 Cs 



Finally, to test the universality of our results, the three- 
body recombination enhancements for the model poten- 
tials of Lennard- Jones type, having either 2 or 3 s-wave 
bound states, are shown in Fig. [6] In this figure, the red 
(black) solid vertical line shows where the <i-wave dimer 
crosses the threshold for LJ with 2 (3) s-waves. The red 
(black) dashed vertical line indicates where the three- 
body state associated with the d-wave level crosses the 
threshold for LJ with 2 (3) s-wave states. The values of 
the s-wave scattering length for such enhancements from 
these two different models differs by only 0.01 r v( jw- This 
small difference is consistent with the small shift of the 
scattering length value where the d-wave dimer crosses 
the three-body threshold (see Fig. [TJ). As discussed be- 
fore, the d-wave scattering length is a better parameter to 
describe d-wave interactions. We have therefore analyzed 
the three-body recombination rate as a function of l/af, 
which is shown in Fig. [7J The agreement between calcu- 
lations using different potential models in fact looks bet- 
ter in this parameterization, where the three-body state 
crosses the threshold at a<i ~ 0.902r v dw for both LJ with 
2 and 3 s-waves bound states. 



This result suggests (although it does not prove rig- 
orously) that the three-body state associated with the 
d-wave dimer is also universal. We also note that at 
the threshold for formation of the three-body state, the 
effective potential relevant for that state does also dis- 
play an inner repulsive barrier [see Fig. 2] (a)] prevent- 
ing atoms from approaching to hyperradii smaller than 
1.5r v dw- This indicates that the properties of such a 
three-body state should also be insensitive to the details 
of nonadditivc short-range three-body interactions. 
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V. SUMMARY 

In summary, we have calculated the three-body recom- 
bination rates for Lennard- Jones potentials with 2 and 
3 s-wave bound states for small positive values of the 
scattering length. Two universal enhancement peaks are 
found at about a — 0.995 r v( jw- hi particular, one of 
the enhancement peaks corresponds to a new universal 
three-body state that is associated with a d-wave dimcr 
and is expected to exist in general for potentials having 
a van der Waals tail. 
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Appendix A: Lennard-Jones potential 



FIG. 8: (Color online) Lennard-Jones potential with A ~ 
0.503r v dw (black thicker curve) and A « 0.414r v dw (red thin- 
ner curve). The horizontal lines show the bound state spec- 
trum with even angular momentum quantum number I. The 
line style solid, dashed, dash-dotted, dash-dot-dotted, and 
short-dashed indicates I = 0, 2, 4, 6, 8 correspondingly. 



The two-body model potentials applied here are the 
Lennard-Jones potential, which can be expressed as a 
function of interatomic distances r: 



V(r) 



r 6 



A^ 



(Al) 



This form of potential has a van der Waals tail —C^/r 6 
at long range and a repulsive core near the origin con- 
trolled by the parameter A. In the present study, the 
values of A are adjusted to give the desired s-wave scat- 
tering length a and number of bound states. Figure [5] 
shows two Lennard-Jones potentials with a zero-energy 
d-wave state and their corresponding spectrum of bound 
state with even I. (The zero energy d-wave state is not 
shown in Fig. [5]). The black thicker curve corresponds 
to A ps 0.503r v dw an d has two s-wave bound states with 
a ps 0.995r v dw- The red thinner curve corresponds to 
A ps 0.414r v dw and has three s-wave bound states with 
a ps 0.983r v dw- It is shown that although the potentials 
at short-range are very different for the two different pa- 
rameters, the scattering length and the spectrum near 
threshold are similar. 



Appendix B: Improved method for calculating 
nonadiabatic coupling matrices 

Equation (|16p is only accurate up to the first order 
in AR. In addition, the value chosen for AR in a re- 
alistic numerical calculation can sometimes be tricky. 
When (R; f2) changes very rapidly, e.g., near a sharp 
avoided crossing, we need to choose a very small step size 
AR. In contrast, when 4> M (R; fi) changes very slowly, 



e.g., at very large distances R, we need to choose a rel- 
atively larger step size AR, or else $ M (R + AR; Q) — 
$ M (R — AR; f2) would be too small, and the accuracy 
would be limited by the machine precision. 

One way to improve the accuracy is to apply the 
Hellmann-Feynman theorem. This theorem can give us 
analytical formulas for the coupling matrices if we know 
the derivative of the adiabatic Hamiltonian ^i? a d from 
the following derivation. First, taking the derivative of 
both sides of Eq. leads to 

[H M (R, O) - U v (R)] (R; fi) = 



$ v (R;Sl). (Bl) 



^(R,n)-±uAR) 



Next, multiplying $ M (R; Q,)* on both sides of Eq. ([Bl 
and integrating over £1 gives 



J dn-f, (R; or [^g ad (R, nj\ <K (R; % 
[u ll (R)-u v (R)] w ' 



where fx ^ v , and 

u v (R) = J dn$ u (R-ny 



d_ 

dR 



OR 



#ad {R, 0) 



®u fan) 

(B3) 

after some manipulation of algebra. The Hellmann- 
Feynman theorem is in principle exact. The matrix ele- 
ments for P 2 can be obtained by 



p 2 - V p p 



(B4) 
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where N c is the number of channels. However, numerical 
studies show that the convergence of with respect 
to number of channels is very slow, making this method 
impractical. 

We now introduce an improved method to calculate 
(i?; f2). It is in principle numerically exact. The 



first hint of the derivation of this method is that Eq. (|B1[) 
seems plausible to be directly solved for -§^^^ (R; fi) by 

^$„(fl;fi)= -[ffadtofiJ-MA)]" 1 * (B5) 
d d 

—H &A (R, ^)-q^U v (R) (R; fi) . 



which can be derived from the normalization condition 
of $, y (i?, f2), as shown in Rcf. [Hj]. Nevertheless, our 
numerical studies show that even without removing the 
singularity, the use of numerical solution packages such 
as "Linear Algebra PACKage" (LAPACK) [H or PAR- 
DISO I13| to solve Eq. (|B6|) directly still gives an ac- 
curate solution x„ (R; f2) in the form of Eq. (|B7j) with 
an unknown c. And once we have the numerical solution 
Xu (R; fi), c can be calculated by 



(B9) 



However, this solution is forbidden since 
Had (.R, 0) - U v (R) is singular: \H ad (R, Q) - U v (R)\ = 
0, meaning that i? a( j (R, fi) — U v (R) is not invertiblc. 
The singularity can also be understood from the fact 
that the equation 



[H ad {R,tt)-U v {R)]xu{R;tt)= (B6) 



does not have a unique solution, Xv {R'i I n fact, any 
functions with the form of 



xAR\V) 



d_ 

dR 



$„ (R; 0) + c$„ 



(B7) 



(where c is an arbitrary number) can be a solution of Eq. 
(|B6[) . The singularity of matrix ff a( j (i?, ft) — U u (R) can 
be removed by considering the additional condition that 



/ 



dCL$„(R,CL)*-^$ u (R,n)=0, (B8) 



Finally, the derivative of <& v (i?; fi) can be written as 



—§ v (R;Si) = Xu{R;to)-c&u(R;n), (Bio) 



which can be inserted into Eq. (fTl"]) and Eq. (fl"3"]) for 
the coupling matrices. The P matrices obtained in this 
way are found to be numerically the same as the one 
calculated from the Hellmann-Feynman theorem up to 
machine precision, proving that our -§ji^y {R; 0) are nu- 
merically accurate. Therefore, this method can give nu- 
merically very accurate coupling matrices. The major 
limitation of this approach is that it becomes computa- 
tionally more expensive to obtain nonadiabatic couplings 
as the number of requested channels increases. Neverthe- 
less, for our typical calculations the number of relevant 
three-body channels is not so large and this method has 
been shown to be quite efficient. 
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